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The permeability of two-dimensional fractures with self-affine fractal roughness is studied via 
analytic arguments and numerical simulations. The limit where the roughness amplitude is small 
compared with average fracture aperture is analyzed by a perturbation method, while in the opposite 
case of narrow aperture, we use heuristic arguments based on lubrication theory. Numerical simu- 
lations, using the lattice Boltzmann method, are used to examine the complete range of aperture 
sizes, and confirm the analytic arguments. 
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I. INTRODUCTION 



The transport of fluids through geological media such as hydrocarbon and water reservoirs involves a combination 
of flow through the microscopic pore space of the rock itself and flows through macroscopic channels such as fractures. 
The first case is relatively well understood, at least in principle by means of models which treat the pore 

space as a random network, and then use effective medium or percolation concepts for the transport. Fractures are 
typically modeled as simple straight-sided channels, with a cubic relation between fluid flux and average aperture, 
and a challenging problem is to understand the dynamics of flow in a macroscopic fracture network |^-|6| . Typically 
the surface of a single fracture appears fairly smooth, aside from some small-scale superficially random roughness, 
and Poiseuille flow in a straight channel is the obvious model for fluid flow. However, more careful analysis |?| shows 
that common geological fractures in fact have correlated, self-affine fractal surfaces. The roughness exponent, whose 
precise definition is recalled below, is usually found [@|| to be close to 0.8, surprisingly insensitive to the material and 
the fracturing process. (Other values may arise from intergranular effects JTof.) 

The aim of this paper is to study the permeability of single self-affine two-dimensional fractures. We shall see that 
it is possible to obtain general analytic expressions in two limiting cases, where the roughness associated with the 
fracture surface is either small or large compared to the mean aperture. The theoretical predictions will be supported 
by numerical simulations, using the lattice Boltzmann method, which of course can address the case of intermediate 
sized roughness as well. In subsequent work, we will consider the fully three dimensional case, and go on to consider 
more general transport processes in fractures, involving both passive tracers and finite-sized suspended particles. 
Complimentary experimental studies on laboratory samples of fractured rock are in progress elsewhere Jll| . 

To the extent that the fluctuations in the height are small compared to the aperture width, the effect on permeability 
is modest, and we will obtain small corrections to the usual cubic law, but the effects on other transport processes are 
much more significant. For example, passive tracer dispersion is very sensitive to the presence of any slow zones in the 
velocity field, and even uncorrelated roughness leads to significant long-time tails |l2|. Similarly, the deposition and 
erosion of solid particles from a surface is sensitive to the local shear stress, which in turn can change dramatically as 
a surface roughens. 

In the opposite extreme of a very narrow fracture, the flow field changes qualitatively as the fluid winds through 
a highly irregular channel. Aside from the challenge of estimating the permeability of difficult geometry in terms of 
its statistical characteristics, there are further effects arising from the fact that the two sides of the fracture originate 
from a single crack. In the fractal case, the spatial correlations between the two sides of a fracture lead to velocity field 
correlations, which again strongly affect tracer motion 13-1^]. The latter references made use of a very approximate 
velocity field, obtained in the lubrication limit, and one of the motivations of this paper is to examine the validity of 
the lubrication approximation for the velocity field. 

The organization of this paper is as follows. We first recall some basic facts about self-affine surfaces, and give the 
algorithm used to generate them numerically, and also explain the lattice Boltzmann method used for the numerical 
simulations. We then consider the case of fractures with a wide average aperture, obtain a perturbative estimate for 
the permeability, and test it numerically. We then turn to narrow fractures, first in the case where the two sides are 
simply displaced normally to the mean fracture plane, and secondly when there is a lateral shift as well. Finally, we 
summarize the results, and indicate the next steps 
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II. PRELIMINARIES 



A. Self-afHne roughness 

We briefly review the mathematical characterization of self-affinity, and its implementation in this paper. We 
consider a rock surface without overhangs, whose height is given by a by a single- valued function z(x,y), where the 
coordinates x and y lie in the mean plane of the fracture. Self-affine surfaces ||l6| display scale invariance with different 
dilation ratios along different spatial directions, remaining unchanged under the rescaling 

x —* \\x y — > A2J/ z — ► A3Z (1) 

Here we consider disordered media, so these scaling laws apply only in an ensemble or spatial average sense. Experi- 
ment indicates that for many materials isotropy can be assumed in the mean plane, implying that there is only one 
non-trivial exponent relating the dilation ratio in the mean plane to the scaling in the perpendicular direction, i.e., 
Ai = A2 = A, and A3 = AS with 

z{x,y) = \~ c z(\x, Ay) (2) 

where £ is the roughness or Hurst exponent |17|| . 

We assume that the process of fracturing the rock is "clean" , in the sense that the rock breaks so as separate along 
one single valued surface without subsequent deformation and without producing loose interstitial material. We shall 
emphasize two complimentary limiting situations, in which the two surfaces are either well separated with respect to 
each other, or very close to each other. If L is the lateral size of the fracture in the mean plane, the typical range of 
the fluctuations (the difference between the maximum and minimum values of z) scales as R ~ L/> . If h is the mean 
width of the fracture, with h <C L, then the two cases correspond to R <C h and R^> h, respectively. 

The local aperture of the crack is the key parameter determining the fracture permeability, and three different cases 
will be addressed. For large aperture, the roughness associated with the two sides of the fracture act independently, 
and it suffices to consider a channel with one rough and one smooth boundary. In the case of a narrow fracture, the 
correlation between the two sides is a key feature, and we consider two possibilities. We first study fractures where 
the upper surface has been simply translated a distance h normal to the mean plane, so that the local aperture a(x, y) 
equals the constant h. Alternatively, the two surfaces may have a relative lateral displacement d in the mean fracture 
plane, accompanied by a displacement h in the perpendicular direction, so that the two surfaces do not overlap. In 
this case the local aperture is given by the random variable. 

a d (x,y) = z(x + d,y) - z{x,y) + h (3) 

It turns out jl8j that d is the lateral correlation length for fluctuations in the aperture, in the sense that a<i(x, y) and 
a,d{x + Ax, y) decorrelate for Ax 3> d. 

In this paper we restrict ourselves to the two dimensional case where the surface is invariant in the y direction, 
z{x,y) = z(x), and the flow is forced in the x direction by a constant pressure drop. In a subsequent paper we 
will extend these calculations to fully three dimensional fractures, but it is convenient, both conceptually and in the 
numerical simulations, to regard the system as having a translationally invariant third dimension. 

B. Numerical Methods 

Our aim is to study various aspects of the fracture permeability which are sensitive to the fracture roughness. The 
first ingredient is the height function z(x, y), a statistically self-affine surface with periodic boundary conditions. The 
periodicity is not a physically essential ingredient here, but has some calculational advantages in alleviating finite-size 
effects. The surface is generated by a Fourier synthesis method, based on power-law filtering of arrays of independent 
random numbers |l{J. The random numbers are generated using a Gaussian distribution, and then modulated by 
an appropriate power law. If Z(k) is the Fourier transform of the initial Gaussian random array, then the Fourier 
transform of the surface elevation is chosen to be 

Z{k) = k-^- x ^Z(h) (4) 

(The 1/2 is appropriate for the self-affine curve used here, and would be 1 for a real rock surface.) 

In most cases the roughness exponent is chosen as the often-observed value £ = 0.8. The amplitude of the roughness 
can be expressed in terms of variance of the height distribution over the full range, 
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a 2 = ((z(x) - (z(x))) 2 ) = I - (z^))) 2 ^, (5) 

which in this case is a 2 = 0.49 ± 0.06 for L = 256. Alternatively, the amplitude is given indirectly in terms of the 
range of the fluctuations as R = RiL 1 * , with Ri — 1.66 ± 0.4. In these equations, and in the rest of the paper, the 
unit of length is that of the spacing between lattice points in the numerical calculations. 

Since we consider fluid flow in a highly irregular geometry, and eventually hope to consider dispersion and the 
evolution of the solid surface due to particle transport, the lattice Boltzmann (LB) method |2(J is particularly conve- 
nient. In this algorithm, fictitious particles move between neighboring sites on a lattice, with suitable collision rules, 
and the boundary of the flow domain is simply a surface of sites where the rule is modified in some way to keep 
the particles out. We use the Face Centered Hypercubic (FCHC-projected) version of the LB model, with a cubic 
lattice in 3 dimensions and 19 velocities (D3Q19 in the terminology of pi]]). The collision operator is approximated 



by a single relaxation time, the Bhatnagar-Gross-Krook model 22 , and the local equilibrium distribution given by 
PH is used. This pseudo-equilibrium distribution locally preserves mass and momentum values, and is formulated 
specifically to recover the Navier-Stokes equation at large length and time scales. For the no-slip solid boundaries, 
we use the simplest implementation of particle exclusion - the bounce-back rule, in which a particle incident on the 
boundary reverses its direction. Periodic boundary conditions are used for the inflow and outflow surfaces. A constant 
pressure gradient forcing the fluid is added in the x direction, while the gap between surfaces extends over z, and the 
geometry is taken to be translation invariant in y. 



III. SMALL SURFACE DEVIATIONS 



In this section we consider the case of a channel with one fractally-rough wall, in the limit where the mean width 
is large compared to the amplitude of the roughness. We (re-)derive an elegant general result for the permeability, 
and then an exact variant which gives the permeability of a perturbed pore space. If the pore space perturbation is 
weak, the leading order correction is easy to evaluate and, when compared to numerical simulation, is seen to provide 
a reasonably accurate estimate of the permeability. 



A. Integral representations 

We begin by deriving an integral representation for the permeability, and then considering a perturbation in the 
boundary shape. Although we suspect that these results are known to many, we are not aware of an earlier publication 
containing them, although the electrical analog of the integral representation is in the literature f23f| . 

Consider a rectangular volume of a porous medium, periodic in all directions, and the following surface integral 



ds ■ 



s 



u(f) p(r) (6) 



Here, u and p are the velocity and pressure fields of a fluid which completely fills the pore space, assumed to satisfy 
the Stokes equations, and the surface S is the complete boundary of the pore space, consisting of the inner grain 
surfaces and the porous regions of the outer boundary. The pressure and velocity fields are periodic, except that the 
pressure jumps by a constant amount P between two opposite faces in one direction, that of the mean flow. Now 
the velocity vanishes on the grain boundaries by the no-slip condition, and there is complete cancelation between the 
opposite faces in the two completely periodic directions, but in the flow direction the remaining two faces of the box 
combine to give 

f kP 2 ^ 
I = P ds-u{r) = -PQ = — (7) 



E 



Here Q is the fluid flux through this end face E, and the minus sign arises because the outward normal is used in /. 
The last equality follows from Darcy's law, where S is the area of end E, k is the permeability, fi the fluid viscosity, 
and L the length of the box. 

On the other hand, applying the divergence theorem to Eq. ^, we obtain: 



I = 



[ dVV-(up) = [ dV {u-Vp)= [i [ dV (u- V 2 u) = -fj, [ dV (Vu) 2 (8) 
Jv Jv Jv Jv 
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where the volume V is just the pore space. We have used first the incompressibility of the fluid, second the Stokes 
equation, and third integration by parts. Comparing the two expressions for /, we have 



(9) 



The analogous formula for the electrical conductivity a, which follows from Ohm's law and the Laplace equation [123] 
is 



0=00 ih! Y ^ ? {di ' 



(10) 



Here <7o is the conductivity of the pore fluid, <j>(r) is the potential and $ its difference across the sample. Note that 
both of these formulae may also be derived by comparing the energy dissipation rate computed microscopically in the 
pore space to the same quantity computed using average fields. 

We next derive an exact variant of the integral expression for the permeability (|J) due to Wilkinson j24j] , which 
allows us to implement a perturbation analysis. Suppose we begin with a pore volume Vq and known Stokes equation 
solutions po and uq, and then contract the volume to V , where the exact (but unknown) solution is p = po + p e and 
u = Hq + u e . The result to follow is true even if the volume change and alterations in the fields are not small, but are 
probably only useful in that limit, hence the suggestive subscript e. Substituting into the dissipation integral on the 
right hand of Eq. [9| and suppressing the summation sign for the moment, we obtain 



dV {d lUj f = / dV (cUo,,) 2 + {diU €tj f 



dV diUQ.j d., 



The last integral on the right hand can be rewritten as 



dV 



di(u j diU Ct j) -Uqj <9j«e,j = / dSi Uqj diU Ct j - dV U J <9jU e ,j 



(11) 



(12) 



Noting that ito, 



-u e j at the surface S due to the no-slip condition, and using fidfu e j = djp € , we may rewrite the 



previous expression as 



dV diU j diU £ 



dV 



(diu € 



(13) 



Using incompressibility, the last term in the integrand can be rewritten V • (up e ) and converted to a surface integral, 
which vanishes if the pressure is held constant on the ends of the porous medium. Finally, we obtain 



L 



S \P 



, , Jv 



(14) 



This result is exact but not useful. However, note that if the relative change in volume is small, O(e), the second 
term is second order in the small parameter, so that to first order the decrease in permeability involves only the 
integral of the unperturbed velocity over the deleted region, which is easy to calculate. 



B. Dependence of permeability on system size 

In this section we will determine the decrease in permeability due to the presence of a self-affine perturbation in the 
lower surface of a straight channel, and in particular its dependence on the size of the channel. To use the previous 
result, we begin with a straight channel of width ho which completely contains the rough-walled one (see Fig. [I]), 
where we have Poiseuille flow 

1 P 

u = — — z(h -z) x (15) 
2fi L 

and an unperturbed permeability ko = h\jYl. The lower boundary is shifted by z(x) > 0, where we assume that that 
range R of z satisfies e = R/h <C 1. The only surviving velocity derivative is 
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1 p 

d z u 0:X = — — (h - 2z) (16) 

so that the dissipation integral is simply 

dV^diUoj) 2 = / dx I dy I dz(d z u , x ) 2 (17) 

Up to first order in e this becomes 

dV{d lU ^) 2 ^ [ — ) Yi^L / -(,-),/,■ } (IS) 




p\ 2 s , r 3 



fij L I h 

where (z) is the average of z over the channel length L, and the cross-sectional area is S = hoY. The mean perturbation 
(z) is computable for a specific profile, but here we wish to use a statistical characterization and relate it to the average 
properties of an ensemble of self-affine surfaces. 

The mean height of the surface is half the range, and therefore has the scaling (z(x)) = ^R\L^ . Substituting, and 
replacing the previous equation in Eq. ^ we obtain the first-order perturbative correction to the permeability 



k sa ko 



(19) 



We have added an adjustable parameter C\, which is expected to lie between 1 < C\ < 2, to take account of the 
distinction between open volume and flowing volume. A value C\ = 1 means that the fracture is equivalent to a 
straight channel of height equal to the mean height of the surface. On the other hand, a value C\ = 2 means that 
the whole region of fluid below the maximum excursion of the surface is not contributing to the permeability and the 
system is equivalent to the maximum channel that does not intersect the surface. 

Note the lowest order correction to the permeability obtained here results from the reduced volume of a fracture 
compared to that of a straight channel enclosing it, and happens to coincide with the lubrication approximation. The 
approach taken here allows us to gain some insight in the drawbacks of this approximation and further improve them. 
Using the previous equation we may write the first order correction to the flow rate, 



Q(L) « Qo 



~ 2ho 



(20) 



To verify this prediction, we consider a pore space consisting of a cubic lattice with periodic boundaries in x and y, 
one impermeable wall at z = ho, and a self-affine rough boundary lying above z = 0. We generate an ensemble of 10 
self-affine surfaces with exponent ( = 0.8, as discussed above. In Fig. |] we display the correction to the flow rate as 
a function of the size of the system L. The straight line is a fit to the expected exponent, which does reasonably well 
except at the smallest and largest values of x. The origins of the discrepancies are first that the discretization used 
in the numerical calculation suppresses any roughness less than one unit, which is a significant fraction of the system 
for small L. Moreover, the system should be periodic in L, but the numerical algorithm used to generate it assumes 
a system size that is a power of 2. As usual, there are computational limitations on the sizes we can simulate, and for 
lengths between 32 and the accessible maximum of 64, we just truncate the periodic surface. Finally, at the largest 
system sizes, we are really outside the range of validity of the estimates. In fact, the second term in the flux in Eq. |2^ 
is greater than 1/3 for the parameters here when L w 16. 

+2.1 

The fitted value C\ given that £ = 0.8 is C\ = 1.30 -1.2, meaning that a portion of the fluid does not contributes 
to convective transport. Expressed in terms of the surface height variance a 1 instead of the span range, the excluded 
fluid represents 50% of the roughness region. The physical consideration absent in the leading order perturbation 
calculation is the fact that fluid velocity sharply decays inside depressions and corners were resistive eddies are likely to 
occur This deficiency is associated with the omitted second term in Eq. |lj, i.e., the correction to the unperturbed 
velocity. In Fig. [I] we magnify two regions close to the lower and upper surfaces, respectively, showing the markedly 
different way the velocity decays towards the surface. 
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We have also simulated a system with an alternative roughness exponent C = 0.95, corresponding to a surface with 
more persistent correlations and in a sense a lower fractal dimension (the dimension of the intersection of the surface 
with a plane normal to it is D = 2 — £ (l6)). In this case the prefactor in the range scaling law is R\ = 0.6 ± 0.12. In 
Fig. H we show the numerical results; again, good agreement with the predicted behavior is obtained. The fit gives 

+3 1 

C\ = 2.04 -i's, suggesting that the low-velocity zones are even more important in this case, presumably because the 
greater degree of correlation enhances them. 



C. Permeability decrease due to zero mean surface fluctuations 



We wish to disentangle the permeability decrease due to volume reduction of the fracture from that due to low- 
velocity zones induced by surface fluctuations. The first step is to consider "zero-mean" self-affine surfaces, whose 
mean height h is constant, but which have a tunable roughness. The previous profile, z{x) for < x < L, is rescaled 
so as to have unit variance, and then multiplied by an adjustable amplitude A, giving adjustable variance A 2 . If we 
assume that the effect of the boundary perturbation is to give low velocity zones which do not contribute to the flow, 
and that the volume of these zones contributing to the dissipation integral is proportional to A, we would write in 
analogy to Eq. EG 



! ZC 2 A 



(21) 



Here A is analogous to R\ and C2 an order-1 fitting coefficient. In Fig. || we show the results obtained by 

numerical simulation corresponding to a system with h = 32 and L = 64 (and, to be specific, Vp = 6.25 x 10 -6 , 
/i = 0.1, Y = 4, which gives Qo = 0.683). The filled circles are the numerical results, and we see the expected 
linear decrease of the permeability with A. In fact, a linear fit to the numerical data gives Qo = 0.72 ± 0.03 and 
C2 = 0.77 ± 0.1, which are reasonable numerical values given our assumptions. 

Further evidence for our arguments concerning the effects of low- velocity zones can be adduced by showing that the 
permeability is unchanged if they are deleted. We begin with the family of adjustable-amplitude self-affine surface just 
discussed, and filter out the smallest wavelengths from the spectrum so as to generate a smother surface. The smooth 
surface is then shifted upwards so as to enclose the deleted surface fluctuations, in the RMS sense. Quantitatively, if 
we have Nl points on the x-axis, with L = iV^A, and filter all wavelengths smaller than A c = ANl/Nc, then 

1 Nl 9 

V u W n=N c +\ n 

The prefactor on the right arises from the discrete form of Parseval's theorem p6fl , and in practice we choose A = 1. 
We then simulate fluid flow in the smoothed fracture, and in Fig. |^ the open squares show the results obtained as a 
function of roughness amplitude A. The agreement between the actual and filtered surfaces is clear. 

The extension of this analysis to the 3-D case is conceptually straightforward, since the decrease in the permeability 
in this limit just corresponds to the decrease in pore space volume. The latter would be L x x L y x R 1 due to the 
self-affine topography of the fracture surface, and the relevant expansion parameter is again R/Iiq. The effects of 
low-velocity regions is presumably the same as well. 



IV. NARROW FRACTURES 



We now turn to the complimentary situation of narrow fractures, where the pore space is winding and convoluted. 
Consider the situation in which a rock is fractured and the two surfaces are simply displaced by a distance h <C L, 
perpendicular to the mean fracture plane, with no lateral shift. The vertical aperture is constant everywhere, but the 
effective local aperture for fluid flow, i.e., the local width of the channel normal to the mean flow direction, strongly 
depends on the local angle between the surface and the mean plane. In Fig. || we show a fracture of length L = 64, 
generated by a self-affine surface of exponent £ = 0.8, separated a distance h = 8, along with the (lattice Boltzmann) 
computed velocity field. 

The theoretical analysis will be based on a kind of local lubrication approximation, where the channel is divided into 
a sequence of quasi-linear segments at varying orientation angles. First, we estimate the length £|| in the direction of 
the mean flow over which the channel formed by the two fracture surfaces can be considered approximately straight, 
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which is a typical size over which the fluctuations in the vertical direction are small compared to the effective aperture 
of the channel. Using the self-affine scaling law for the correlation function, 



2C 

(23) 



oi(^) = {[z(x,y)-z{x + Z h y)] !> ) = 
where £ is a microscopic length, say a grain size, and (f>(£) is then on the order of a grain size squared, we estimate 

z(x,y)-z(x + ^,y)r.t^j (24) 

A segment of length ^ is roughly straight when this quantity is a small fraction of h, which yields 



i/C 



Returning to the entire fracture, each £|| -channel is oriented at some angle 9i with respect to the mean plane, 
and has effective aperture cij = hcosdi, and length £| = £||/ cos^. Assuming Poiseuille flow in each ^-channel, the 

corresponding local permeability is fcj = af/12. We compute the pressure drop across the fracture by adding the the 
local pressure drops along the straight £||-channels along the whole path. Noting that the flow rate Q — Uh — Uidi 
is constant, 



N N 



i=l i = l 

N 



£ 

; = 1 



12/i 



h 2 cos 2 8i J \h cos 0i J \ cos 9 



= — £ cos (**) = — £ cos 

Finally, noting that N — L/^ 3> 1 is the number of channels, we can convert the sum into an average over the 
distribution of angles, and write for the permeability 

h 2 1 

k = 12 <cos-4(0)) (2?) 

We can give a simple if heuristic estimate of the permeability as follows. Since the exponent £ < 1> the channels 
have small vertical fluctuations, and we can approximate the cosine as 

cosfl= , ^l-l(^M) 2 (28) 



Substituting in Eq. |?], we obtain 
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(29) 



A more convincing evaluation of the angular average makes use of the actual height distribution. Experimental 
measurements indicate a Gaussian distribution for the spatial correlation in heights, Z = z(x,y) — z(x + £m) [ jl8| , 
and in fact our numerical procedure for generating self-affine surfaces also gives a Gaussian distribution for Z . We 
illustrate this point in Fig. ^, where we plot our generated probability distribution function of Z, corresponding to 
different values of £||, along with their Gaussian fits 
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We also plot the (numerically-obtained) dependence of the mean spatial correlation, (Z 2 ) 
choices of £, confirming the scaling given above in Eq. |2^. 
The angular average is then given by 



(cos" 4 9) = 



1 



dZ p{Z) 
g*(gj|) 

£11 



CT z(£||)> on £|| for two 



(31) 



Using Eqs. |23|,^, this result is consistent with the previous estimate Eq. ^9| based on simple scaling arguments. 

Moreover, aside from the purely numerical coefficients, we can argue that scale invariance requires the angular 
average to have the form given in (j3l|). The self-afhne property of the fracture geometry implies that the probability 
distribution of heights p(Z) should only depend on the rescaled variable Z/^f, Jl8fl, More specifically, 



p(Z) = a;\a)^[Z/a z (a)] 



(32) 



where the prefactor comes from the normalization. If (cos" 4 6) is evaluated for this distribution, we obtain a variant 
of Eq. where the numerical coefficients 2 and 3 are replaced by the moments of the function tp, leading to the same 
scaling form for the permeability. 

Finally then, using the leading term in Eq. [Jl] for the angular average, and Eq. ^5] for the value of we have the 
result for the permeability of a narrow two dimensional self-affine fracture, 
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2l. 2 

c 



(33) 



The principal approximation used in obtaining this result is that the fracture aperture may be regarded as a sequence 
of almost-straight segments. 

To test Eq. |3^ numerically, we first recast it in terms of the fluid flux. For a straight channel of height h, width Y 
and length L, and pressure drop P, we have flux Q = h 3 PY '/12p.L. Thus 



PY 2-2C 5C-2 

Qo-Q~C 3 — i— h— 

D/i 



(34) 



where C3 is another adjustable parameter expected to be of order one, and we have taken i to be the (unit) lattice 
spacing. 

In Fig. ^| we present the correction to Qo as a function of the distance between the opposite surfaces h, for both 
roughness exponents C = 0.80 and ( = 0.95. We find the good agreement for the predicted exponent for the case 
C = 0.80. The only adjustable parameter is the coefficient, which is found to be C3 = 3.1 ± 1.2 in fairly good agreement 
with the expected value. The fitted exponent is 2.43 ± 0.08, which agrees with the predicted value (5£ — 2)/£ = 2.5. 
On the other hand, in the case C = 0.95 the numerical results deviates from the asymptotic behavior at small distances 
between the surfaces. That is because the successive powers given by the expansion of Eq. 33 are very close to each 
other (the difference between them is (2£ — 2)/£ m 0.1). Again if only the coefficient is fitted, we find C3 = 0.35 ±0.21. 
Fitting also the exponent, but using only the numerical simulations with large separation between surfaces, we obtain 
2.71 ± 0.39 where (5C - 2)/C = 2.89. 

In Fig. [?] we compare numerical results corresponding to systems of different sizes with the same roughness exponent 
C = 0.80. It can be seen that the accuracy in the exponent is improved as the size of the system is larger. We obtain 
a value 2.52 ± 0.04 for the largest system with size L — 256. 

As mentioned, the crossover between small and large surface roughness, (compared to the mean aperture of the 

5C-2 

fracture), can be addressed numerically. The narrow fracture regime, where Qo — Q cx h 5 , is expected to be valid 
when £|| <C L. The critical mean aperture value h c which corresponds to £|| ~ L turns out to be, using the fitted value 
of C3, h c ~ 100. On the other hand, the small surface-deviation scaling behavior can be obtained from Eq. (|2f~ 
where L is the constant size of the system and ho is now the variable mean aperture of the fracture h, 



Qo-Q 



PY 3Cii?iL c , 2 
h 

Yip, 2 



(35) 



In Fig. U we display the correction to Qo in the range 2 < h < 128 for a system size L — 64. A small deviation 
from the narrow fractures regime can be observed, starting at h = 64. As the size of the gap becomes larger, the 
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correction to the unperturbed flow rate Qq decreases as a power of h. Therefore, both larger computational time 
(starting the simulation from an initially at rest fluid the time interval to reach the steady asymptotic value grows as 
the square of the gap size p7j ) and higher accuracy in the fluid field close to the solid surface are needed. Our current 
computational limitations prevents us to simulate systems where L 3> £11. Nevertheless, we believe that the observed 
deviation correspond to a crossover towards a smaller exponent (h 2 ) as expected from the previous discussion. 

V. FRACTURES WITH SHIFTED SURFACES 

In the field, when a rock is fractured its two sides may be shifted laterally parallel to the mean plane by geological 
processes. We now consider how our results for permeability are modified by this effect. The situation is somewhat 
different in two dimensions than in three. In order that a two dimensional fracture remain open it is necessary that 
the two sides do not touch, whereas in three dimensions in the presence of a lateral shift it is easy for the two surfaces 
to touch at one or a few points in the plane, while the aperture is still open to flow. Furthermore, in two dimensions 
when the fracture nearly touches at a point, the permeability will be dominated by the large pressure drop needed to 
force fluid through this narrow gap. In fact, in this case an estimate of the permeability is simply k w a^ lin /I2, where 
a m i n is the minimum value of the aperture. In this section we address the complimentary case where the fracture is 
distinctly open. 

As discussed in Section 2, in the presence of a lateral shift d, the fracture aperture is now the random function 
ad(x, y) given by Eq. || rather than the constant h. However, there is a residual correlation between the two sides of 
the fracture, which allows us to relate the permeability to the self-affine statistics. First, we determine a condition 
on the shift for the fracture to be open. Provided the average properties of the surface are statistically stationary 
(independent of a;), we have (od(x)) = h, and the variance of the aperture distribution is 

a 2 a (d) ee (M*) - M.x)>] 2 > = ([z(x + d) - z{x)f) = al{d) = cf>{i) (|V (36) 

The aperture is surely open when a a (d) <C h, which gives d <C ((h/i) 1 ^, but from Eq. ^5|, this is just d <C £m. In 
other words, the shift must be small compared to the typical length of a straight segment of unshifted channel, and 
there is little change in the geometry compared to the unshifted case. We expect, therefore, that the previous result 
Eq. HH should apply at least up to a different value of the numerical coefficient C3 . 

In Fig. ^ we show LB simulation results for the correction to the unperturbed flow rate Qo, for different values of 
the shift d. We see that the behavior is consistent with these arguments for all values of the shift, with deviations 
occurring when the shift is too large for the mean aperture. Several other shift values were simulated obtaining the 
same general behaviors as in those shown in Fig. ^. As is often the case with this type of scaling argument, reasoning 
in terms of asymptotic limits leads to conditions that one quantity must be much larger than another, but in practice 
the range of validity is much wider. 

VI. CONCLUSIONS 

We have studied the permeability of two-dimensional self-affine fractures, using asymptotic analytic arguments for 
wide and narrow apertures. Numerical simulations using the lattice Boltzmann method have verified the predictions, 
and also suggest a smooth crossover between the limits used in the derivations. We have obtained expressions for the 
permeability in which the usual expression for straight channels is modified by terms related to the Hurst exponent 
characterizing the fracture surface. 

In the wide-gap case, we obtained a perturbative correction related to the roughness exponent and amplitude. 
The result is formally identical to the lubrication theory prediction, but the corrections are known, in principle, and 
furthermore we can understand the discrepancies between the derived numerical coefficient and simulation in terms 
of the effects of low-velocity zones which do not contribute to transport. In the narrow gap case, we also derive an 
analytic correction related to the roughness exponent and amplitude (or microscopic length £) . The analysis is based 
in a local lubrication approximation, dividing the whole fracture into a chain of approximately linear channels with 
varying orientation angles respect to the main flow direction. The results is equivalent to a correction due to the 
tortuosity of the fracture, where the tortuosity itself depends on the gap size. Finally, we show that in the case when 
there is a lateral shift between the surfaces, there is no qualitative change compared to the unshifted case. This 
is basically due to the two dimensional geometry, where in order to have an open channel the shift must be small 
compared to the size of the typical length of a quasi-linear segment of unshifted fracture. 







The two-dimensional nature of the problem imposes restrictions on the range of validity of analytic estimates, but 
the numerical results are in general agreement over a much wider range. 

The extension of this work to fully three-dimensional fractures is in progress. The case of a wide gap does not require 
any further significant conceptual effort. However, the narrow gap case the fracture geometry is significantly more 
complicated because it is not feasible to use any quasi-one-dimensional approximation for the flow paths. Numerical 
simulations are certainly feasible, at least up to size limits imposed by one's computational hardware, but additional 
ideas are required for analytic arguments. Future work will consider diffusion, and also convection plus diffusion, in 
these self-affine fractures, in both two and three dimensions. The LB method is quite simple for numerical simulations, 
and we will explore the way in which the roughness exponent enters into the quantitative behavior, along the lines 
given here. 
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FIG. 1. Example of the geometry and velocity field in a fracture with one self-affine surface of roughness exponent £ = 0.80. 
The enlargements show the difference in the velocity decay near smooth and rough boundaries. 

FIG. 2. Perturbation in the flow rate in a wide, rough channel as a function of the length L of the system for £ = 0.80 and 
C = 0.95. 

FIG. 3. Flow rate as a function of the relative amplitude of the surface roughness. Circles and squares correspond to the 
fully rough and the small- wavelength filtered surfaces, respectively. 

FIG. 4. Flow field in a narrow trough fracture with a constant gap, and exponent £ = 0.80. The enlargements illustrate the 
effect of the effective aperture on the flow field. 

FIG. 5. Top: semi-log plot of the distribution of heights Z — z(x) — z(x + £u) for different values of £m for a self-affine surface 
with exponent £ = 0.80. The solid lines are a Gaussian fit to the numerical results. Bottom: variation of the mean spatial 
correlation (Z 2 ) with £||, for ( = 0.80 and < = 0.95. 

FIG. 6. Change in flow rate as a function of gap h between vertically-displaced self-affine surfaces, for exponents £ = 0.80 
and C = 0.95 and length L = 128. 



FIG. 7. Change in the flow rate as a function of gap h for £ = 0.80 and various L. The solid line is a fit to the largest 
(L — 256) system. 

FIG. 8. Change in the flow correction as a function of the gap h for £ = 0.80 and L = 64. The solid line is the to the narrow 
fractures regime. At large gaps (h > 64) a small deviation towards a smaller exponent can be observed 

FIG. 9. Change in flow rate correction as a function of gap h for exponent ( — 0.80 and various lateral shifts d. The solid line 
is a fit to the case d — shown in Figure ^. The two vertical dashed lines divide the regions (Jill) < (d/f) c and (h/t) > {d/£) c , 
for d = 22 and d = 46. 
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